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Abstract 

We present a Schwinger-boson approach for the RVB state of the spin-1/2 Heisenberg 
antiferromagnet on a triangular lattice. It is shown that a Gutzwiller projection of the 
mean-field state that includes both antiferromagnetic and ferromagnetic decouplings leads 
to optimizing the RVB pair amplitudes within a self-consistent approximation. The re- 
sulting state yields, by Monte Carlo simulations, energies and spin-spin correlations in 
excellent agreement with the exact diagonalization result on finite lattices (up to 36 sites). 
We conclude that the optimized RVB wavefunction possesses a long-range three-sublattice 
order. 
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The spin-1/2 antiferromagnetic Heisenberg model on a triangular lattice has been a 
subject of great interest in recent years. Unlike the square lattice case, the correspond- 
ing Ising model is disordered at all temperature with finite entropy This was also the 
original model for Anderson's proposed resonant valence bond (RVB) state 0. The sys- 
tem was thereafter often refered to be in quantum liquid state for the lack of long-range 
magnetic order (LRMO). However, this simple picture was not so obvious: The quantum 
fiuctuations may as well lift the classical degeneracies and lead to ordered ground-state. 
Many studies were further inspired by the magnetic properties of doped Cu02 planes in 
the copper-oxide superconductors where the frustration again plays a central role. Al- 
most every possible scheme has been employed in the literature. These include the earlier 
variational approaches[0, 0, |], the spinwave and large- theories[^ 0, the spontaneous 
symmetry breaking consideration [0], the numerical series expansion methods]^, H, and the 



exact diagonalizations on finite clusters |[10|, [TT], |T^. Most of recent analytic aproaches tend 
to conclude with LRMO. But they all seem to suffer from uncontrolled approximations or 
inaccuracies. On the other hand, the numerical groups are divided into with and without 
LRMO (or just right on the critical point as suggested by Singh and Husep]), mainly due 
to insufficient knowledge on larger clusters and different methods of extrapolation. The 
main question, whether this frustrated system possesses LRMO remains controversial. 

Even RVB states could have LMRO if the bonds decayed sufficiently slowly. This was 
first found on a square lattice by Liang, Doucot and Anderson|T^. Does this also take place 
on the the triangular lattice? In this Letter, we shall address to this question based on a 



new use of the Schwinger-boson formalism|T|, 15|. The central idea is to work with the 
RVB states resulting from Gutzwiller projection of mean-field states onto the right Hilbert 
space. This method has been working extremely well on the square latticefl^. Note only 
the ground-state energy and the staggered magnetization obtained are virtually identical 
to the exact ones (cited by|T^), the low lying spin-flip spectrum also agrees with the best 
numerical result found by supercomputer (cf. |]17|). On the triangular lattice, one needs 
to face the frustration, we shall show that it is important to have, at the mean-field level, 
both the antiferromagnetic and ferromagnetic fields in order to account for this frustration. 
With a proper balance between the two, the Gutzwiller projection then leads to optimizing 



2 



the RVB state. The latter is also explicitly derived within a self-consistent approximation 
in the "loop-gas" statistics [l^, |13|. Our approach provides a simple but powerful way of 
founding the best RVB state and the right mean-field solution regardless of the frustration. 
The optimized parameter-free RVB state is found to have energies much closer to the exact 
ones than those previously found [0, ^. The system is predicted to be long-range ordered 
with a staggered sublattice magnetization = 0.275. 

The nearest neighbor antiferromagnetic Heisenberg model may be written in terms of 
Schwinger bosons as|T^ 



H 



E 

<ij> 



s 



(1) 



(2) 



where 

Aij = bi^bji - biibj^, Fij = ^ bUja, E ^^t^ta = 2S. 

We have deliberately broken up in (|l|) the Hamiltonian into two pieces (the first two terms): 
The first (second) is suitable for an antiferromagnetic (ferromagnetic) mean-field decou- 
pling. The separation {u is a free parameter to be fixed) is motivated by the fact that 
the pairing of nearest neighbor spins is frustrated on the triangle lattice. A ferromagnetic 
field will therefore be essential for an accurate description of the original system. This will 
become clear below when we minimize the energy of the RVB state. The partition function 
of the system can be calculated via 

Z(/5) = Tr[PGP] = Tr|[nA]p|' 



i=l 



where p = exp{—(3H), and Vg (Pi) stands for the Gutzwiller projection operator for the 
whole lattice (the ith site) which enforces the constraint in (Q). Introduce the so-called 
coherent state, \ai >= exp (^X](7=Ti ^L^io-) |0 >, where = (6j|,6j|) are ordinary complex 
numbers. We can further write [ITli 



■ N 

.1=1 



< {«I}lpl{«i} > 



(25)! 



E 



_d d_ 

db~dbL 



{at 
2S 



=0} 



(3) 



(4) 



provided that the matrix elements of p are known. 

Eq. d^) would be rigorous if the exact matrix elements were used. This is of course 
rarely possible for interacting systems. One therefore needs a trial p. We shall consider in 
this Letter a mean-field type of approximation in which the Hamiltonian is replaced by 

^mf = ^0 + A ^ hlK + E + QiAj] + h.c] . 

i,cr <ij> 

The Lagrangian multiplier A and the mean fields Dij and Fij are to be determined via the 
usual self-consistency requirements, 



J2 < blAa >= 1, 



n 



l + iy 



Q^ 



1 - Z/ 



< > . (5) 



mf ^ft^r Fourier transforms can be readily diagonalized by a Bogoliubov transformation. 
This results in a mean-field spectrum ci;k + Mk with 



CUk 



Mk = (Q_k - Qk)/2, 



Ak = A + (Q_k + Qk)/2. 



The matrix elements of p can be calculated accordingly. We find|14 



< {"I}|Pmfl{"»} >= ^o(/3) exp ^ 



where 



iy^') = 5kexp(-/?Mk) 



W^f = (DkM)5kSinh(/3c^k) 



5k = [(Ak/cc^k) sinh(/5a;k) + cosh(/5cUk)]-^ (6) 

Note that Qk = <5k> = -^-k- 

Now return to (^) and consider in what follows 5 = 1/2. Write Z{P) = Zq{I3) x YJv 
and take a closer look at the structure of 1^. The differentiations amount to dividing the 
whole lattice into subsets of self-avoiding loopsfl^, |15[. One then sums over all possible 



configurations of the loops. This kind of loop gas statistics has been known[13, 18| for the 



RVB states. The more general case|IJ, is that the bond, say, from sites i to j is replaced 
by a 4 X 4 transfer matrix Tjj of the form (cXz below is the third Pauli matrix) 



/ 



G 







OzGija^ 



(7) 



It brings, upon taking Pi, see (||), the prefactor (xa&iT + Xi2&l| + Xisbii + XiM^) to a new 
one at site j via multiplying T^- to the column of the x's- ^//^ respectively, the 

(inverse) Fourier transforms of Tracing over the matrix (divided by two) after 

completing a loop gives its contribution. 

To illustrate this prescription, one can decompose Yat at a given Jq into (all j'^'s below 
are self-avoiding) 

N-l 



^^=E E V}v-„-i({j4) X [Tr(G,„,, 



G 



(8) 



The arguments of Y^^n-i represent the sites excluded. Various physical quantities can be 
calculated by multiplying the operators to p and then compute the matrix elements as in 
(I). In the coherent states one has the simple substitutions, b]^ bi^ {d/db\^), 

which may modify the transfer matrix at site i. In general, operators that conserve boson 
numbers can be presented by some 4x4 matrices inserted at the relevant sites. One then 
can proceed to evaluate the modified as usual. For the spin-spin correlation Sj ■ Sj we 
have the matrices 





This leads us to {i ^ j) 



1 N~l ( on 

< ■ S, >= — E 5: YN-n-i{{3k}) X - E Tr[a.G51a,G 



(2) 



G 



3m~l3m 5 



/~<(2) _ /~( ^^ 



(9) 
(10) 



We shall take here a specific mean-field solution Consider the simplest symmetry Qij = 
Q, Dij = iD exp(i3(/>jj) with being the polar angle of the bond. Let 7k = (1/3) cos(k- 
e^), Fk = (l/3)E-=iSin(k- e,) where ei = (1/2,-^3/2), = (1/2, +^3/2), and eg = 
(— 1, 0). At zero temperature, the mean-field dispersion cUk/A = cuk = y^(l + (ii7k)^ — {d2Tky 
becomes gapless. Namely we have (2; = 6 is the coordination number) 



1 + (ii7±ko = ±c?2r±ko, 



di = zQ/X, 



d2 = zD/X. 



(11) 
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where ko = 27r/3(l, -\/3). Condensation, denoted by a which has to be inserted in 
occurs at ±ko [or their equivalences on the reciprocal lattice spanned by Gi = 7r(2, -\/3) 
and G2 = (0, 47r/v^)]. Consequently, the mean-field equations of (H) read ("BZ" means 
the first Brillouin zone) 



a 



1 r VSd^k 

2 Jbz 



1 + c?i7k 



D 



l + v 

Q_ 

V 



1 



OL 1 

2^'^" + 4 7i.z 



OL V 

2 + 4 iBZ 



87r2 
87r2 



7k(l + (ii7k) 

CJk 



(12) 



The matrix elements P^j/s can be found upon substituting the the solution back to (^, see 
(0) below. 

We next show that this set of VT^'s minimizes the RVB ground-state energy at z/ = 
within a self-consistent approximation for the loop gas. Going back to the transfer 
matrix (|^), w'^^ vanishes at zero temperature. The physical picture of (0)-(|Tnp become 
apparent: Here stands for the amplitude of a RVB pair connecting site i to j, more 
explicitly, iyij(6||(i]| — [the superscript "(2)" is hereafter dropped]. Fat is nothing 

but < '^'kVbI'^RVB ^ with |$pyB ^ containing all possible superpositions of the pair 
products. Eq. (P) then gives the explicit way of computing the spin-spin correlations. The 
question is reduced to finding the set of best amplitudes that minimizes the energy. 

Rigorous evaluation of the RVB expectation values appears problematic. We thus in- 



troduce below a self-consistency approach. Let us first approximate jY^^n-i 



V 



n+l 



This assigns a uniform weight for a loop of (n + l) sites. We are still left to deal 

with the self-avoiding restriction. Let us try to ignore this in the first place, denoting the 
corresponding y by y^. The correct y is then obtained by properly identifying the over- 
counting. This allows a full analytic summation over the loops using Fourier transforms. 
The matrix after summing over the paths connecting i to j reads [i 7^ j, see (|10|)], 

/ 



R,; 



1 exp(— ik ■ Fjj) 
- |W^k/2/o|' 



(13) 



1 VTkVl/o 
— 1) = 1 [cf. (|15D below]. Now y can be 



The parameter y^ can be determined by Tr(Rj, 
approximately recovered by grouping the extra winding at a given site (say 1) into Xjy 



which gives 



y-^ = y^' X TtRu/2 = (3/2) x y^ 



-1 



(14) 



More detail on this over-counting problem has been dicussed in [113 



The same procedure can be applied to the spin-spin correlations. We now have one path 
starting from i to j and the other from j to i. Ingoring again the self-avoiding restriction 
and the overlaps between the two paths, the result of the summation is simply Rjj and Rjj. 
However, taking two independent paths results in an additional over-counting of a factor 
3/2 [see (p^j which has to be deducted. Sutbstituting (|1^) back to (§) and assuming, as in 
the mean-field case, the three-fold symmetry for W^, we end up with the expression (drop 
yo hereafter as is variational) 



E- 



bond 



4" 1 - 



7k 



1 >p T^W^ 



subject to Wk = = ~W-k and the constraint 

N ^1- 



3 
2' 



(15) 



Eq. (|T5D can be treated via the method of Lagrangian multiplier. For simplicity, one can 
assume here real W^s. It is now a straightforward exercise to obtain (picking up the right 
solution satisfying \Wk\ < 1) 

1^2 rk 



with 



1 + rfi7k + V(l + rfi7k)2 - (d2rk)2 



di (l/iV)Ek[7k/(l-W^k') 



(16) 



(17) 



d2 (l/iV)Ek[rkW^k/(l-W^k')]' 
Note that ([T6|) is nothing but the zero-temperature limit of (^) when the mean-field solution 

discussed above is used. Moreover, it can be easily checked that (0) and ([17|), after some 

manipulations using the specific form of (|T^, reduce to (|TT]) and ( [T^ ) of the mean-field 

result with u = 0. Note that, in passing to the continous limit, the same condensation a 

is also required in (^) and ([l^)- This proves our claim that two approaches give the same 

RVB state. 
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lim < Si ■ Sj >- 

\rij\^oo 



Solving (|T|) and or equivalently (|nD and (0), we find a = 0.275, di = -0.7223, 
d2 = 1.572 and -E'bond ~ —0.1899. Because of the condensation at the particular modes 
±ko, the system possesses a long-range three-sublattice order with 

i, 7 G same sublattices 

(18) 

-a /2 i,j G different sublattices 
The structure of (|^) and (^) allows also a direct evaluation via Monte Carlo simulation. 
This is important for justifying our theory. It amounts to sampling over the distribu- 
tion of the loops. At zero temperature, the rule to calculate the spin-spin correlations is 
particularly simple: For a given configuration, 

—3/4 < ipili^R > i,j& same loop n{i) — n{j) = odd 
+3/4 < ipili'R > hi ^ same loop n{i) — n{j) = even 
hj^ different loops 

where n{k) stands for the order of the site k in its own loop. Each loop contains even number 
of sites. In the simulation, the loop configurations are updated by randomly choosing a 
pair of nearest or next nearest neighbor sites and exchanging their loop connections with 



< ipilSi ■ SjiijR >-- 



a probability satisfying detailed balance condition ||13|, p!5|] . For the triangular lattice, the 
sign problem shows up in that some of the loops contribute negative weights. We use their 
absolute values to regulate this problem. This limits the practical sizes of the simulation. 
The periodic boundary condition for the RVB amplitudes also presents a subtle problem. 
Normally, we should use those obtained by summing over a discrete set of k's in accordance 
with the finite size of the lattice. If we did so, there would be no sign problem for sizes up 
to 6 X 6. But the energy would also get higher. We therefore use the amplitudes calculated 
on a much larger lattice (200 x 200). We then go to the smaller size and cut them into a full 
hexgon (i.e. with inversion symmetry). Bonds outside are moved into the hexgon via the 
periodic boundary condition except for those that are required to vanish by symmetry [e.g. 
for i = (0, 0) and j = ±(n, —n) when the two axes are chosen to be (1,0) and (—1, -\/3)/2]. 
The correlations as a function of distance can be obtained by averaging over the six-fold 
equivalences. 

We wish to compare our results simulated on finite lattices with those based on exact 
diagonazations. On a 4 x 4 lattice, the energy per bond is -Ej^yg = —0.1754(3) vs. 
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-^exact ~ —0.1782 (cited from |Tl|). The difference is probably due to that the symmetry 



of the bonds are not very suited on this cluster. The result on a 6 x 6 lattice is presented 
in Table 1. One sees that they are in excellent agreement. Perhaps, the only notable 
difference between the two sets is at = 3. The result also shows a nice approach to (|T8D , 
i.e., a three-sublattice magnetization. We stress that our wavefunction is parameter-free. 
These results themselves should be of interest for variational approaches in the presence of 
frustration. Finally, the small difference in the approximate and the true values of -f'bond 
should justify our analytic work. 

Let us summarize the main results of this Letter. Starting with a Schwinger-boson 
mean-field state that includes both the antiferromagnetic and ferromagnetic fields, we have 
shown that a Gutzwiller projection leads to optimizing the RVB state on the triangle lattice. 
This in turn justifies the use of the two fields. The ground-state energies of the optimized 
state were calculated rigorously on finite lattices and were found to approach virtually 
the exact values. It is apparently much better than all previous vaiational states|P, ^. 
we have demonstrated, though our prediction for larger lattices is still based on analytic 
approximation, that by properly optimizing the RVB pair amplitudes one is most likely to 
get long-range correlations. Namely, the system is predicted to possess a long-range order 
within the RVB context. 

The author is indebted to J. L. Lebowitz, Z. Y. Weng and T. K. Lee for helpful sugges- 
tions. He also thanks the hospitality of Rutgers University. This work was supported by 
the National Education Committee of China, and by the US Air force Office of Scientific 
Research through Grant No. 92-J-0115. 
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Table Caption 

Table 1 The spin-spin correlation C{n) as a function of distance on the 6x6 lattice. The 
simulation takes up to 2 x 10^ runs and the statistical errors are about ±0.002. The 
exact result is quoted from Leung and Runge[|T^]. 
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n = 


1 


2 


3 


4 


5 


6 


RVB 


-0.186 


0.155 


-0.062 


-0.065 


0.118 


0.117 


exact 


-0.1868 


0.1535 


-0.0548 


-0.0664 


0.1136 


0.1174 



Table 1 



